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The optical and electronic characteristics of devices based on GaAs 
(leds, laser diodes, etc.) are adversely affected by the dislocations 
originating in the substrates. We demonstrate by means of thermoe- 
lastic analysis that the primary cause for the observed dislocation 
density patterns in Czochralskipulled GaAs single crystals, which 
serve as a source for substrates, is crystallographic glide, induced by 
the excessive thermal stresses arising during the growth process. 
First, we formulate a tractable model for crystal growth. We obtain 
the temperature distribution in the crystal by solving the quasi- 
steady-state partial differential equation for heat conduction subject 
to appropriate boundary conditions. The closed-form solution in- 
cludes time, pull rate, axial location, radius, convective and radiative 
heat transfer coefficients (h r + he), and a fixed ambient temperature 
(T a ) among the variables. Next, from the temperature profiles we 
determine the radial, tangential, and axial stress components acting 
on the GaAs boule. These stresses permit the evaluation of the 12 
resolved shear stress components for the {111}, (110) slip system. We 
postulate the sum of the absolute values of the 12 components (o tl „) to 
be proportional to the dislocation density within an additive constant. 
Employing a to i as a parameter, we have constructed dislocation 
distribution contour maps for {100} GaAs wafers which are in good 
accord with the dislocation patterns observed on KOH-etched wafers 
cut from near the top end ofCr and Te-doped GaAs boules. A detailed 
examination of the effect of the numerous parameters on the dislo- 
cation density of Czochralski-pulled GaAs is also given. Only by a 
drastic increase of T a and a substantial decrease of h r + h c would 
one be able to overcome the natural limitations imposed by the 
thermal and mechanical properties on dislocation density. Finally, 
we pay attention to the effects of elastic anisotropy and interfacial 
heat flux, discuss the philosophical and mathematical difficulties 
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associated with finding a true transient solution, and provide some 
practical suggestions. 

I. INTRODUCTION 

A number of recent papers have shown that the performance of 
GaAs-based devices is adversely affected by dislocations. For example, 
Brantley and Harrison 1 have observed that the degradation rate in 
diffused electroluminescent diodes increased by an order of magnitude, 
and was accompanied by dislocation generation, during forward bias 
aging when a compressive load was also applied. Subsequently, Zaes- 
chner 2 has found that external stress alone is sufficient to cause diode 
degradation and gave a quantitative account of the change in light 
output with time in terms of the kinetic properties (multiplication and 
velocity) of dislocations, considered to be nonradiative recombination 

centers. 

Even in the absence of deliberately imposed forces, a large enough 
density of grown-in dislocations in GaAs alters device behavior. In 
particular, very recently Roedel et al 8 have correlated the reduction in 
external quantum efficiency of Si-doped GaAlAs leds with increasing 
dislocation density and established that these defects act as nonradia- 
tive recombination centers. Moreover, the dislocation density in the 
GaAlAs epitaxial layer essentially duplicated that of the GaAs sub- 
strate. Therefore, in view of its importance, we have undertaken an 
investigation of the primary cause for the generation of dislocations in 
GaAs substrates grown by the Czochralski technique and then attempt 
to employ this knowledge in suggesting growth conditions which 
facilitate the elimination or at least reduction of these defects. 

Some plausible mechanisms by means of which dislocations can be 
incorporated into growing GaAs boules are (i) dislocation propagation 
and multiplication from an imperfect seed, (w) vacancy condensation 
into dislocation loops, and {Hi) crystallographic glide relieving exces- 
sive thermal stresses. It has been recognized relatively early in the 
course of semiconductor crystal growth that thermal stresses may lead 
to slip and dislocation generation. In 1955, Billig 4 discovered that the 
etch pit density of Ge wafers, obtained from pulled ingots, increased 
with the magnitude of the imposed temperature gradient. Further- 
more, the pits (representing the dislocations) were distributed in a 
definite pattern and slip bands were also observed. On the basis of 
these experiments, Billig concluded that high thermal stresses* gave 
rise to slipping and dislocation generation. He also made an order-of- 
magnitude estimate of the thermal stress to offer a more quantitative 



* The Stress was given by the product of the thermal expansion coefficient, Young's 
modulus, and the radial temperature drop. 
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underpinning to the postulated mechanism. The required temperature 
difference between the core and the edge of the crystal was obtained 
from a standard heat transfer calculation for a semi-infinite cylinder in 
the steady state. The base of the cylinder is at the melting point, T f , 
while the sides are dissipating the heat by convection. 

Almost simultaneously with Billig, Bennett and Sawyers 5 of Bell 
Laboratories found a hexagonal "star pattern" of pits along definite 
lines on etched Ge slices cut from (111) pulled ingots and gave the 
following qualitative description of the thermal stress effect: The heat 
enters the growing Ge crystal at the solid-liquid interface and leaves 
through the other surfaces by radiation and convection. Hence, each 
cross section of the crystal must have a cooler periphery than core 
and, consequently, on account of thermal contraction, the periphery 
must be in tension and the core in compression. If the resolved shear 
stress components of the {111}, (110) slip system resulting from this 
tension are sufficiently high, then the resulting plastic deformation 
would exhibit the symmetry of the observed dislocation pattern. 

In 1958, Penning 6 reported that etch pit arrays in germanium can 
also be introduced by radial heat flow during slow cooling from 850° C. 
The patterns show a crystallographic orientation-dependent symmetry 
and the dislocation density is highest at the edge, intermediate in the 
center, and lowest along an internal annulus. Following the earlier 
work, Penning also interpreted his results in terms of thermal stresses 
and formulated a semiquantitative model. In his view, it is reasonable 
to assume that two alternative avenues of stress relief are open to the 
crystal. In one case, the strain is entirely plastic and is completely 
relieved by the generation and motion of dislocations. In the other 
case, the thermal stress is mostly elastic, but a small constant fraction 
is released by plastic flow; hence, the dislocation density corresponding 
to one of the 12 {111}, (110) slip systems is proportional to the amount 
of slip, governed by the appropriate shear strain component. A com- 
parison of the observed and estimated etch patterns favors the path of 
incomplete stress relief. 

Unfortunately, Penning 6 withheld the details of his assumptions and 
treatment for a promised forthcoming paper. To the best of our 
knowledge, that paper has not been published in the intervening 20 
years. However, it seems clear from the brief account given in Ref. 6 
that his calculations employed a simple parabolic temperature distri- 
bution in the crystal and the predicted etch figures were based on the 
shear components of the thermal stresses. In our view, even today, 
Penning's conceptual framework presents the most fruitful starting 
point in crystal-growth related thermal stress analysis. 

A strong confirmation of plastic deformation taking place in semi- 
conductors during the growth cycle was provided by the work of 
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Lederhandler, 7 who correlated infrared birefringence with etch pit 
density in <111) Si. This investigation shows that the source of 
birefringence is the elastic-plastic deformation occurring under the 
influence of a radial temperature gradient. The frozen-in stress seen 
by birefringence is acquired by the crystal near room temperature and 
is equal and opposite to the thermal stress relieved by dislocations and 
plastic flow during growth. According to a more recent study, similar 
to germanium, the dislocation density in silicon exhibits a minimum at 
a location intermediate between the core and the lateral surface. 8 

Although the dislocation density in III-V compounds is occasionally 
orders of magnitude greater than in germanium and silicon, there is a 
paucity of information on the source of grown-in dislocations in the 
semiconducting compounds. In the case of liquid-encapsulation Czo- 
chralski-(LEC) grown GaP, Nygren 9 has found that, consistent with the 
thermal stress mechanism of dislocation generation, their density 
increases toward the periphery of the ingot; he also observed traces of 
slip. Moreover, the "frozen-in" stress was determined on the transpar- 
ent crystal by a photoelastic technique showing equal tangential and 
radial stresses at the core and a disappearance of the radial stress at 
the periphery. 

However, there is a lack of consensus with respect to the genesis of 
dislocations in GaAs. Based on their results using a modified Grem- 
melmaier-type magnetic puller, Steinemann and Zimmerli 10 claimed 
that thermal strain is inconsequential in causing dislocation generation. 
But they admitted that the omission of a heat shield from their 
apparatus made the growth of 1-cm diameter, low-dislocation density 
GaAs impossible even when a long neck was employed to eliminate 
the dislocations arising from the seed. Under these unfavorable con- 
ditions, the dislocation density increases toward the external surface. 

Subsequently, Brice 1112 has identified the partial pressure of arsenic 
(As 4 ) as an important contributing factor to the observed dislocation 
densities in GaAs. Crystals grown by the horizontal Bridgman 
method," and in a syringe puller 12 at pressures between 0.8 and 1.3 
atm, show a corresponding monotonic increase in their dislocation 
densities, probably as a result of gallium-vacancy condensation into 
dislocation loops. 11 Nonetheless, Brice has also suggested 12 that the 
radial temperature gradient via the thermal stress mechanism is re- 
sponsible for the enhanced dislocation density near the edge of his 
crystals and surmised that, at lower pressures, the same mechanism 
may be responsible for the entire distribution. This latter notion is 
supported by Plaskett et al 13 who examined, by etching and x-ray 
topography, the slip line patterns of GaAs crystals grown by the 
horizontal Bridgman method at pressures somewhat below 1 atm. 
They concluded that the dislocations were created by plastic defor- 
mation arising from thermal stresses. 
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Therefore, in view of two decades of accumulated experience in the 
area of the melt growth of semiconducting crystals, it is reasonable to 
propose that the initial formation of dislocations in GaAs crystal pulled 
from stoichiometric melts by the Czochralski technique is primarily 
due to the excessive thermal stresses associated with the growth 
process. A major objective of this paper is to test this proposition by 
the formulation and analysis of a tractable model for growth which 
describes the dislocation distribution in terms of material and growth 
system parameters, thus permitting a direct comparison with experi- 
mentally derived dislocation patterns. As a first step, we require 
realistic temperature profiles for a growing ingot. A simple steady- 
state heat transfer model introduced by Brice 14 has given the experi- 
mental temperature distribution in germanium 14 and ZnW0 4 , 15 ob- 
tained by a thermocouple embedding technique, in a consistent man- 
ner. He assumed that the growing crystal can be represented by a 
stationary cylinder, the base of which is held at T f , while its lateral 
surface and top dissipate the heat into a medium at a constant 
temperature by convection and found the solution in the classic 
monograph of Carslaw and Jaeger."' However, this solution is only 
valid for a fixed ambient temperature and excludes growth rate as a 
parameter. In the present work, we have been able to include growth 
rate in the solution because a moving boundary quasi-steady-state, 
partial differential equation for heat conduction is considered instead 
of a steady-state one. In a future paper, we shall report the exact 
solution of the steady-state partial differential equation subject to 
convective heat transfer from the sides of the cylinder into an ambient 
with a linear temperature profile which abruptly changes slope at some 
height. Obviously, that result closely corresponds to the thermal ge- 
ometry of lec growth wherein a sharp break in temperature occurs at 
the B 2 0.j-gaseous ambient interface. 

The quasi-steady-state temperature profiles are used in deducing 
the radial, tangential, and axial thermal stress components for a 
growing cylinder in a closed form. Then these stresses are employed in 
evaluating the 12 resolved shear stress components of the {111}, < 110> 
slip system causing glide. Invoking Penning's hypothesis, 6 the dislo- 
cation density is taken to be proportional to the sum of these shear 
stresses. Computer simulation of the theoretical results is facilitated 
by some numerical techniques. Among these, we concentrate on the 
inversion of the formula for shear stress as a function of radius and 
angle into a polar plot of dislocation density contour lines. 

Next, model calculations for the dislocation density patterns are 
presented, which show the effect of a realistic variation in the heat 
transfer coefficient, radius, pull rate, and time. Furthermore, the 
theoretical contours are compared to the dislocation density patterns 
exhibited by KOH-etched wafers. To ascertain that the GaAs crystal 
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encountered plastic flow above the critical resolved shear stress, con- 
sideration is given to the actual magnitude of the resolved shear stress 
in terms of the best estimates of the physical and geometrical param- 
eters appropriate in the growth system. In addition, we discuss the 
possible effect of elastic anisotropy, the magnitude of the axial heat 
flux, and the philosophical difficulty of finding a true transient as 
opposed to a quasi-steady-state solution. Finally, some important 
conclusions are enumerated and practical suggestions are given to aid 
in the lowering of the dislocation density in GaAs. 

II. THEORY 

2. 1 Quasi-steady-state partial differential equation for heat conduction 
during Czochralskl growth 

In general, severe mathematical problems are encountered in per- 
forming heat transfer calculations involving a change of state. 1 - 1 
However, by an idealization of the Czochralski growth process, we can 
formulate a realistic but tractable model which permits the relatively 
uncomplicated determination of the time- and growth-rate-dependent 
temperature profiles prevailing during crystal growth. The following 
set of simplifying assumptions are introduced: 

(i) Initially, at t < 0, the semi-infinite space between z = and -« 
is completely filled by a stoichiometric liquid solution of gallium and 
arsenic (melt) contained in a crucible and held at the melting point of 
GaAs, T f (see Fig. 1). 



z 



s = o- 




-z = o 



--Z = -pt, 



CRUCIBLE 
CONTAINING 
Ga/AsMELT 



Z = -pt 2 




GaAs 
CRYSTAL 




-S = 



Fig. 1— Stationary (z) and moving (s) coordinate systems for Czochralski growth. The 
stationary and moving systems are anchored at the top (seed) -end and the solid-liquid 
interface, respectively. Above the interface, the temperature of the surroundings is 
uniformly T„. 
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(ii) The crystal is grown by the removal of the crucible in the -z 
direction at rate p, and its shape is cylindrical from top (seed) to 
bottom (tail). 

(Hi) The growing boule is surrounded by an ambient fluid at tem- 
perature T a = constant < 7>, while the Ga/As melt is continuously 
maintained at Tf. 

(iv) The planarity of the solid-liquid interface at T f is essentially 
unchanged by the dissipation of the heat of fusion of GaAs. 

(v) As the top of the crystal always remains stationary at z = and 
the ^olid-liquid interface moves to z = -pt, the sense of growth is 
opposite that of the practical situation. This choice of direction is 
dictated by mathematical convenience, and it leaves the heat transfer 
unaffected as long as T a is a constant. 

(vi) The temperature distribution is not significantly influenced by 
crucible or crystal rotation. 

(vii) At the top and on the lateral surface of the boule, the heat loss 
(flux) is proportional to the temperature difference between the surface 
and the ambient fluid. This is otherwise known as convection boundary 
condition governed by Newton's law of cooling. 

Figure 1 sketches the coordinate system relevant to our Czochralski 
growth model. In the interest of clarity, the vertical axes (z and s) have 
been displaced from the crystal axes. In a cylindrical coordinate system 
(r, 0, z; t) with the origin at the center of the stationary crystal top, the 
partial differential equation for heat conduction takes the form 

dT (d 2 T IdT d 2 T 



where »c(cm 2 /s) is the thermal diffusivity. The diffusion equation is 
independent of 6 on account of the cylindrical symmetry of the 
boundary conditions. 

It is convenient to transform eq. (1) into a coordinate system (r, 9, 
s, t) embedded in and moving with the solid-liquid interface to facili- 
tate the application of the boundary condition T = 2} at that location. 
This can be readily accomplished 18 by observing (Fig. 1) that 

s = z+pt (2a) 

r = t. (2b) 

Noting that, from eqs. (2), ds/dz = 1, ds/dt = p, dr/dz = 0, and dr/dt 
= 1, the rules of partial differentiation yield 18 

d 2 T_d 2 T dT _ dT dT 

dz 2 ds 2 Ht'P'ds ~dr 

Substituting eqs. (3) into eq. (1), we have the partial differential 
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equation with respect to moving coordinate axes 

d 2 T ldT d 2 T pdT ldT ... 

— 5- + + — r---=- + -TT- W 

dr r dr ds k ds k Br 

According to Rosenthal's theory for moving heat sources, 19 an observer 
moving with the coordinate system fails to detect any change in 
temperature with time in his surroundings. Accordingly, dT/dr = and 
eq. (4) reduces to 

d 2 T ldT B 2 T pdT ._, 

dr 2 r dr ds z k ds 

The preceding equation is called the partial differential equation for 
the quasi-steady state (qss) 1819 and is expected to be valid after a 
crystal of some length has grown. The solution of the boundary value 
problem can be written in a simpler form by introducing the following 
dimensionless parameters: 

p. I, f_£, #-i. (6) 

r r r 



Then eq. (5) becomes 

d 2 T ldT d 2 T_ dT 

where 

Pi =pr /K. 

We attempt to separate the variables by substituting the product 
function 

T = ffWRipWW) (8) 

into eq. (7). The exponential factor in eq. (3) is based on Smolu- 
chowski's solution for field-enhanced diffusion 20 wherein the differen- 
tial equation is formally similar to eq. (7). Using eq. (8), eq. (7) 
separates into 

d 2 R 1 dR p\ 1 d 2 * 



Rdp 2 + P Rdp 4 * cty 2 



(9) 



where a is the separation constant. 

It can be readily seen that the R differential equation, describing 
the radical variation in temperature, is given by 

(fl,{ dH +i* = 0, (10) 



d(ap) otp dap 
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which is Bessel's differential equation of order zero. 21 Bessel functions 
of the first kind and order zero, Jo(ap) satisfy eq. (10) so that 

Rip) = Map). (11) 

The differential equation for the axial temperature distribution is 

/? 2 *--^7F = 0, (12a) 

where 

/? 2 = « 2 +^. (12b) 

Obviously, the exponential function is a solution of eq. (12a). Hence, 
we can write, in general, 

* = A sinh fi(b -M+B cosh P(+ t - yfr), (13a) 

where 

r radius 

and A and B are constants. 

Since the value of the separation constant a is so far unrestricted, 
the complete solution of the qss problem is obtained by means of eqs. 
(8), (11), and (13) in the form of the infinite sum 

T = T a + e"*' 2 £ MpoLn)[A n wdh Pntyt - i>) 

n-\ 

+ BnCO&k Pnty, - $)]. (14) 

2.2 Boundary conditions 

According to the previously outlined Czochralski growth model, at 
the top of the growing ingot and along the lateral surface, the heat is 
dissipated by Newtonian cooling. Expressed mathematically, this 
means that, at the cylindrical boundary, 16 

dT 

— + h(T-T a )\ r=ro = (15a) 
of 

or 

dT 

— + hi(T- r a ) p -i = 0, (15b) 
dp 

where h = hi/r is the heat transfer coefficient. On the stationary top 
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surface of the cylinder, 



or 



— + h(T - T a )\s= pl = (16a) 

ds 



^- + /n(T-r a )|^, = 0. (16b) 

By 

Applying eq. (15b) to eq. (14) leads to 

^ (a nP ) + h Ma n p) | p =, = 0. (17) 

dp 

An important property of the Bessel functions of the first kind is the 
recursion formula 16,20 

£^ =-</,(*), (18) 

ax 

where Ji is the Bessel function of order one. Rewriting eq. (17) in 
terms of eq. (18) provides the n characteristic equations 

-OLnMctn) + hMa n ) = (19) 

for the eigenvalues a„. A combination of eqs. (16b) and (14) results in 
a relationship between A„ and B„ of the form 



Bn 
fin 

where 

Pi 



A n = h p ^, (20a) 



hp^ + h,. (20b) 

Substituting eqs. (20a) into eq. (14), we obtain for the temperature 
profile of the growing ingot 

T = T a + e p ^ /2 i Jo(a nP ) ^ [/* P sinh £„(* " ^ 

+ P„COSh finty, - Ml (21) 

where the a n s are the eigenvalues of eq. (19). 

We determine the remaining constant B n from the boundary condi- 
tion 

T = T f (22) 
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at the planar solid-liquid interface (\p - 0). Then eq. (21) reduces to 

GO D 

Tf-T a = £ Jo(a nP ) -f- [/ipsinh M, + /?„cosh M t ]. (23) 

n=l Hn 

To find 5„, we must expand the constant Tf — T a as a series in Jo. 
Based on the orthogonality properties of Bessel functions, subject to 
the characteristic eq. (19), it can be shown that an arbitrary function 
f(p) can be expanded by means of the Bessel series 16 



f(p) = I KManp), (24a) 



where 



2a 2 f 1 

Kn = <h*J- 2v"r2, \ Pf(p)Jo(<XnP) dp. (24b) 

(hi + an)do{a„) J 

Putting f(p) = Tf — T a and in view of another recursive property of 
Bessel functions, 16,20 

*£*> -*/„<*), (25) 

ax 

eq. (24b) integrates to 

K _ 2a n JMn){T ( - T a ) _ ZhATf- T a ) 

(h] + an)Jo(dCn) (hi + Ctn)Jo(a n ) ' 

Term-by-term evaluation of B n is accomplished by a combination of 
eqs. (23), (24a), and (26). Then we find 

B m 2hAT f -T a )(3 n x 1 (27) 

(/if + al)Jo(a n ) Apsinh /?„»//, + /?„cosh fijfa * 

Finally, substituting eq. (27) into eq. (21), the qss temperature 
distribution for Czochralski growth becomes 

T ~ Ta = 2h ie "^i ManP) 

yn\ t- a n )do\a n ) 

hn sinh R„(\b, — \b) + R-c.nsh RJ,b. — \L) 

(28) 



Tf - T n „=, (hi + a'n)Jo{a n ) 

hpsinh finii'i - ty) + /?„cosh /?„(>//, - ip) 



/i,,sinh /?„i//, + /?„cosh /?„»/>, 

where the summation is over the n eigenvalues of eq. (19). It is 
convenient to repeat here the previously given definitions of the 
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variables and parameters occurring in eq. (28). These are 

, , Pro . Pi . 

/ii = /ir , /?i = — , ftp - — + Ai 

/ft. = ' 



t-(a+f) 



p-i.f.i + g and *-2. (29) 

2.3 Thermal stress 

The temperature distribution induces a thermal stress field in the 
growing cylindrical ingot as a result of spatially inhomogeneous ther- 
mal contraction. The appropriate stress components can be obtained 
from classical thermoelastic theory. 22 In essence, in this theory an 
additivity hypothesis of elastic and thermal strains is superimposed on 
Hooke's Law. There are only few exact solutions of the thermoelastic 
equations. However, for a long isotropic cylinder with axisymmetrical 
temperature distribution, an exact description of the stress components 
is possible if the displacement is only radial, the center suffers no 
displacement, and the lateral surface is free of traction. 22 This set of 
assumptions provides the so-called "plane strain" solution which is 
subsequently adjusted by means of Saint-Venant's principle to take 
into account the absence of traction at the top and bottom ends. 
Consequently, the following final expressions have been derived for 
the radial, o r , tangential, o g , stress components: 

— r^Gif**-?! Iydr| (30a) 



•-^(af B -* + ?i>*- J, i (30b) 



aE (2 



«-—,{nL Trdr ~ T y (30c) 

where E, v, and a are Young's modulus, Poisson's ratio, and the linear 
thermal expansion coefficient, respectively. 22 

To calculate the components explicitly, it is convenient to rewrite 
eqs. (30) in terms of the dimensionless variable p and integrate with 
respect to a n p. Then we have 
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aE 1 



Ta n pd(a n p) 5 Ta n pd(a n p) 



(31a) 



aE 1 



o e = 



l-val 



Ta n pd{a n p) + -5 
P 



Ta n pd(anp) - alT 



(31b) 



a z = 



a£ 1 
1-pSI 



Ta n pd(a n p) - a\T 



(31c) 



All the integrals in eqs. (31) can be readily performed by applying eq. 
(25) to eq. (28). Hence, the stress components become 



a r = 2h\e 



j>i*/a 



aE 



I- v 



(T f -T a ) 



hpsinh finitt - $) + /3/iCosh Pntyt - *P) 



„ti (/i? + al) Jo(a n )(h p sinh fi^f t + j6„cosh Mt) 

J\{a n ) Ji(otnp) 



OtnP 



(32a) 



a B = 2hie 



.Pxt/Z 



aE 

1- v 



(T f -T a ) 



" /ipsinh fintyt - ^) + ft„cosh /?„(»//, - ^) 
A (/i? + a 2 n )J (a n )(hpSinh M, + /? n cosh M,) 



X 



J\(OL n ) Jl(<Xnp) _ 

+ Jo(a n p) 

Oin Otnp 



(32b) 



a z = 2h 1 e^ /2 ^-(T f -T a ) 
1 — v 

y hpsinh /U^ - jj + /? n cosh ff w (fr - jj 
„.i (/i? + a^)Jo(a n )(/i P sinh /3„i//, + /?„cosh /?„»//,) 



2Ji(a n ) 



- e/o(a R p) 



. (32c) 



It should be noted that eqs. (31) are not changed by an additive 
constant {T a ). 
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2.4 Resolved shear stress 

The major mechanism by means of which dislocations are first 
introduced during the Czochralski growth of GaAs is crystallographic 
glide caused by the excessive thermal stress. According to Schmid's 
Law, 2:i glide occurs when the resolved shear stress {ors) exceeds a 
certain critical value, the so-called critical resolved shear stress. Since 
the dislocation density, d,, is proportional to the glide strain, we may 
also take it to be proportional to the o RS within an additive constant. 

It has been previously shown that slip in GaAs is associated with 
the {111}, (110) slip system. 24 This is also the case for other semicon- 
ductors crystallizing in the diamond or zinc blende structure. 25 The 
{111}, (llO) slip system represents 12 permissible glide operations, 
that is, four {111} slip planes, each containing three possible <110> slip 
directions. In this section, we give o RS for each of the 12 distinct slip 
systems, obtained from the principal thermal stresses o r , og, and o z . 
The calculation is shown in detail for one system; the remaining 11 
can be derived by analogous procedures and only the final results will 
be quoted. 

The coordinate system appropriate to perform the stress transfor- 
mations is presented in Fig. 2. To begin, it is necessary to determine 
the stress components acting on the xy, xz, and yz coordinate planes. 




Fig. 2— Coordinate system for thermal stress transformations. For a crystal growing 
in the [001] orientation, the direction of the radial (o r ), tangential (o fl ), and axial (o 2 ) 
components of the thermal stresses are shown. One of the 12 resolved shear stresses 
(o X! ) is illustrated, which acts on the (111) slip plane in the [110] slip direction. 
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By standard tensor transformation, 20 we find 

o x = Or cos 2 6 + o g sm 2 (33a) 

o y = a r sin 2 + o e cos 2 (33b) 

o xy = (o r — o )sin cos d (33c) 

a* = a*. (33d) 

As an illustrative example of the procedure followed, we calculate the 
ors acting on the (111) plane in the [110] slip direction. The coordinate 
system of the ors is the x', y', z' system wherein z' and x' are parallel 
with the [111] normal and [110], respectively. Then ors becomes 

ors = Oxz- = a*cos[.r';t]cos[z':*:] + a y cos[.x';>']cos[2 / ;)'] 

+ o z cos[x' z]cos[z' z] + o X y(cos[x'x]cos[z'y] 

+ cos[x' y]cos[z' x]) , (34) 

where the extended notation has been used and the bracket signifies 
the angle between the indicated axes. In the [001] crystal growth 
direction, it is easy to evaluate the direction cosines; they are sum- 
marized as follows: 

cos[ ] x y z 



x' -\/2/2 V2/2 

y' 

z' 73/3 V3/3 V3/3 

A combination of eqs. (33), (34) and the table provides 

(111), [110] ors=- — (o r - o g )(cos 2 d- sin 2 0). (35) 

fa 

In a like manner, all the o RS can be evaluated. Introducing the 
abbreviation 

o r = Or- og (36a) 

and 

a z = o g - o g (36b) 

and taking advantage of some standard trigonometric identities, the 
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following complete set of ors is obtained: 



Slip 
Plane 


Slip 
Direction 






ff«S 




(111) 


[110] 


a r cos 20 

6 


(111) 


[oil] 


V5 

6 


o z -a,- sin 2 6 - -£ sin 20 




(111) 


[ioi] 


x/6 
6 


o 3 -<7 r cos 2 - -^ sin 20 


(in) 


[ioi] 


6 


<j r cos 2 -a 2 - -y sin 20 




(in) 


[oil] 


76 

6 


o z -a r sin 2 6 + ■£ sin 20 




(in) 


[110] 


a r cos 20 

6 


(iii) 


[OH] 


76 
6 


o z -a r sin 2 - -^ sin 20 




(iii) 


[ioi] 


76 
6 


a r cos 2 -Oz + -^ sin 20 




(iii) 


[lio] 


76 

CTrCOS 26 

6 


(ill) 


[OH] 


76 
6 


a 2 -a r sin 2 + -^ sin 20 




(ill) 


[101] 


76 
6 


a r cos 2 -a a - -^ sin 20 

Zt 




(Hi) 


[iio] 


V 
1 


6 

— 


tx r cos 20. 





(37) 



Inspection of eqs. (37) shows that only five of the 12 stress equations 
are independent. Therefore, we define five ors functions differentiated 
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by Roman numeral subscripts I through V and transform eqs. (37) into 



Slip 


Slip 


Plane Direction 


(111) 


[110] 


(111) 


[HO] 


(111) 


[lib] 


(111) 


[110] 


(111) 


[Oil] 


(111) 


[on] 


(111) 


[Oil] 


(111) 


[OH] 


(111) 


[101] 


(111) 


[101] 


(111) 


[101] 


(111) 


[101] 



ORS 



^6- . 

oi = — — a r cos 20 

6 



ou = 



OlU = 



6 
6 



2 

o z —o r —7= sin 6 sin (8 + tt/4) 
V2 



o z —Or — == sin 6 sin (0 — it/ A) 

V2 



aiv = — 



ov = - 



6 

6 



o 2 —Or —= cos 6 sin (6 + 7r/4) 
v2 

_ _ 2 

o z +o r — = cos sin (0 — 7r/4) 

V2 



(38) 



where we have also employed trigonometric angle-sum relations. 

Finally, since slip occurs regardless of the sign of ors, we take the 
dislocation density to be proportional within an additive constant to 
the quantity a tol which is defined as the sum of the absolute values of 
the 12ors's acting in the {111}, <ll0> slip system. Such correlation 
clearly implies that the stresses are mostly elastic which are not 
completely relieved by plastic flow. 6 Therefore, we have the propor- 
tionality 



gL oc atot = 4 | oi | + 2[| on | + | am | + | orv | + | av |]. 



(39) 



Of course, a„„ can be explicitly obtained by sequential substitutions 
from eqs. (38), (36), and (32) into eq. (39). 

* 

III. COMPUTING METHODS 

To investigate the effect of material and growth system parameters 
on the dislocation distribution in GaAs ingots grown by the Czochralski 
technique, the closed form solutions for the temperature profile [eq. 
(28)], thermal stress components [eqs. (32)], and the sum of the 
resolved shear stresses [eq. (39)] must be evaluated and displayed in 
suitable graphical forms. An essential function in these equations is 
the Bessel function of the first kind orders zero (Jo) and one (Ji). In 
principle, they can be obtained from Taylor series representations; 21 
however, for large arguments the convergence of the series is very 
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slow. Fortunately, Abramowitz and Stegun 27 provide useful algorithms 
for both Jo (their eqs. 9.4.1 and 9.4.3) and Ji (their eqs. 9.4.4 and 9.4.6) 
with an accuracy of better than 1(T 7 . 

Summation of the Bessel functions over the complete spectra of the 
eigenvalues a n is required by all our results of interest. The character- 
istic equation [eq. (19)] immediately suggests not only that a„ is a 
function of hi but also that if hi -* 0, a n is the root of Ji{a n ) = 0, while 
as hi -► oo, a n is the root of J (a n ) = 0. Therefore, the characteristic 
equation has as many roots as the Bessel functions themselves, i.e., 
their number increases to infinity. Since the first root ai at hi —> 
occurs at Ji(«i) = 0, ai — > 0, it is necessary to derive the limiting 
expression for ai for small hi from a Taylor series expansion of eq. 
(19). Then we find 

which is very accurate up to hi =* 0.1. For larger values of hi and a, up 
to i = 6, a numerical table for the roots of the characteristic equation 
is given by Carslaw and Jaeger. 16 These results can be adapted for 
computer use by parabolic least-square fits in hi* which we find to be 
accurate within better than 100 parts per million up to hi < 3 . 

In view of the inverse ai dependence observed in our key expres- 
sions, typically six terms were judged to be adequate to obtain conver- 
gent sums in a„. If, however, functional values very near the solid- 
liquid interface are of interest, a certain caution is warranted. It can be 
readily shown with reference to eq. (28) that, if \p t is large (long boule), 
the radial part of the function is multiplied by the approximate 
quantity e~ M . At the top of the crystal \f> is sizable; hence, as /?„ 
increases with n, e~ Pn * readily tends to zero. In contrast, if v// is small 
(near the solid-liquid interface), many additional e~ Pn * terms are 
needed to achieve convergence. 

The preceding considerations permit the evaluation of the temper- 
ature profiles in the growing ingot according to eq. (28). Since 7>is a 
constant value while T a is an adjustable parameter, it is more conven- 
ient to present the complement of T — T a /Tf — T a , i.e., 

T f -T a T f -T a 

so that the temperature is always referred to a true constant. In Fig. 
3, Tf — T/Tf— T a is shown as a function of the dimensionless radius 
p = r/ro and distance -z/r Q = \p, - \p from the top of the crystal at h 
= 0.3 and 0.6 cm -1 . It is advantageous to employ -z/r Q for the axial 



* A sixth power for a, and a second power fit for a» through a B were used. 
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0.2 0.3 0.4 0.5 0.6 0.7 0.8 

RELATIVE RADIAL DISTANCE, p 



0.9 1.0 



Fig. 3— Relative temperature drop (T, - T/T f - T„) versus relative radial distance 
(p = r/r ) for GaAs at h = 0.3 and 0.6 cm ' [eq. (28)]. The labels are the distances in 
centimeters from the top of the crystal. Other relevant parameters are t = 3600 s, p = 
0.001 cm/s, length = 3.6 cm, r = 2 cm, and k = 0.04 cm 2 /s. Multiplication by T, - T a 
provides the temperature decrease with respect to the melting point (1238°C) of GaAs. 



coordinate since it designates at all times an invariant location in the 
ingot. It can be seen in Fig. 3 that increasing h leads to enhanced 
cooling at the same position and also to a sharper radial temperature 
gradient toward the periphery of the crystal. 

The radial, tangential, and axial thermal stresses, obtained from eqs. 
(32), normalized to (aE/1 - v)[(T f - T a )/200] are presented in Fig. 4, 
as functions of p and -z/r at h = 0.6. In accord with the qualitative 
observations in Section I, the stress components are compressive at 
the core, while the outer surface is in tension (a g = o z > 0, 07 = 0). The 
stress components in Fig. 4 give rise to a to t the sum of the absolute 
values of a RS acting in the {111}, (110) slip system, which is taken to 
be proportional to the dislocation density of a (100) crystal. In Fig. 5, 
we show the radial and axial variation of a t()/ normalized to (aE/1 — 
v)[(T f - T a )/200] at h = 0.6 in the (100) as well as (110) directions. It 
should be noted that in consonance with experimental results for Si 8 
and Ge, 6 there is a density minimum at p ^ 0.6 and the core is more 
perfect than the periphery. Moreover, apart from the core, the dislo- 
cation density is always higher along the (100) than the (110) direc- 
tion. 

Although Fig. 5 provides an acceptable graphical representation of 
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0.1 



0.3 0.4 0.5 0.6 0.7 

RELATIVE RADIAL DISTANCE, p 



0.8 



Fig. 4— Radial (o r ), tangential (o»), and axial (a*) stress levels versus relative radial 
distance (p = r/r ) for GaAs at h = 0.6 cm ' [eqs. (32)]. Multiplication of the levels 
by (aE/1 - v)[(T f — T„)/200] yields the stresses in absolute units. The labels are the 
distances in centimeters from the top of the crystal. Other relevant parameters are t = 
3600 s, p = 0.001 cm/s, length = 3.6 cm, r = 2 cm, and k = 0.04 cm 2 /s. 

the thermal stress effect on dislocation distribution in some directions, 
it does not allow a ready visualization of the distribution over an entire 
wafer surface. Obviously, an explicit plot based on eq. (39) can either 
illustrate atot as a function of p at a constant angle (as in Fig. 5) or as 
a function of at a series of fixed ps.* What one thus needs to serve as 
a ready standard of comparison with experimental information on 
etched wafers is a polar plot of p versus at constant levels of a to t. 
Since an exact inversion of eq. (39) in this manner is impossible, a 
numerical formulation to give p as a function of atot at a constant 
must be attempted. The inversion method adopted takes into account 
the fact that any curve in Fig. 5 possesses a slowly descending lower 
branch up to p a 0.6 and a more rapidly increasing upper branch to p 
= 1. Along the upper branch p = /"(a^; ft 0, h,p,t= constant) can be 
represented by a parabolic least-square fit of degree four. The excel- 



* 6 = and 45 degrees correspond to the (100) and (110) directions, respectively. 
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PARABOLIC FIT 




RELATIVE RADIAL DISTANCE, p 

Fig. 5— Sum of the absolute values of the 12 resolved shear stress levels (o, ot ) or 
dislocation density versus relative radial distance (p = r/r ). The variation in o (ot is 
shown in the (100) and (110) directions of an <001> wafer. Multiplication of the levels 
by (aE/l - v)[(T/ - r„)/200] yields the stress in absolute units. Relevant parameters 
are h = 0.6 cm-, t = 3600 s, p = 0.001 cm/s, length = 3.6 cm, r = 2 cm, and k = 0.04 
cmys. 



lence of the fit is shown by the superimposed dots along the explicitly 
evaluated curves in Fig. 5. The same type of numerical procedure is 
unsatisfactory along the lower branch due to the very slow change in 
CTtot with p up to p ^ 0.2. Therefore, in the lower segment, the p value 
corresponding to a set a„„ was obtained by linear interpolation within 
p intervals of 0.025. 

In practice, p versus a t „t was evaluated at up to 20 equally spaced 
stress levels in 2.5-degree intervals. By interconnecting the points 
corresponding to the same constant a,,,,, a polar plot (p versus 6) of the 
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stress or dislocation density contours was constructed. Without dwell- 
ing in detail on the cumbersome logical choices demanded by the 
development of the plotting routine, at the very least, it should be 
pointed out that when occasionally near the minimum (p « 0.6), at a 
given atot, the parabolic fit to the upper branch yielded a p value 
smaller than the extrapolated one from the lower branch, the point 
corresponding to the lower branch was utilized. 

In Fig. 6, we present the polar plot of the dislocation distribution for 
the top wafer of a (100) GaAs crystal pulled by the Czochralski 
technique. The growth conditions and the temperature and stress 
profiles are identical with those for -z/r = and h = 0.6 in Figs. 3, 4, 
and 5. The contour lines shown are equally spaced and normalized to 
(aE/1 - v)[{Tf — T„)/200]. As an artifact of the stress profile compu- 
tation in = 2.5-degree intervals, the upper and lower branches of the 
closed contours near the center of a quadrant do not intersect in a 
point as they theoretically should, but are interconnected by brief 
straight-line segments. 

It can be seen that the dislocation distribution exhibits fourfold 
symmetry and that the <110) (6 = 45 degrees) is a mirror direction. In 




TOP, 3.6 cm 
LONG 
3600 sec 



<100> 



Fig. 6— Constant o lul or dislocation density contour lines for the top wafer of a <001> 
GaAs boule at h = 0.6 cm -1 . The otot levels are labeled, which can be converted into 
absolute stresses when multiplied by (aE/1 - v)[(T, - T„)/200]. The appropriate 
parameters are t = 3600 s, p = 0.001 cm/s, length = 3.6 cm, r = 2 cm, and (C = 0.04 
cm 2 /s. 
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all directions, there is a density minimum at p = 0.6. Moreover, an 
absolute minimum in the pattern prevails along (110). In moving from 
the (110) to the (100) direction (6 = 0), an increase in the density is 
observed. The crowding of the contours near the periphery of the 
crystal implies a sharp gradient in dislocation density, in accord with 
the previous experimental work summarized in the Introduction. 



IV. RESULTS AND DISCUSSION 

4. 1 Estimation of the input parameters 

Inspecting eqs. (28), (29), and (32), one can immediately realize that 
CTtot and the resulting dislocation density pattern are affected by a 
number of parameters and independent variables. Consequently, it is 
necessary to obtain reasonable estimates for the values of these input 
parameters, in advance of a detailed examination with respect to their 
effect on o tot . In Table I, the range and/or value of these sundry 
quantities of interest are summarized, except for the heat transfer 
coefficient. Fortunately, the thermal and elastic properties are rela- 
tively insensitive to temperature change in the range near T f . The data 
sources for Table I and the methods of evaluation and extrapolation 
will be given elsewhere. 28 

Perhaps the most important parameter is h, which is the constant of 
proportionality for Newtonian cooling on the cylindrical boundary and 
at the top of the crystal [eqs. (15a) and (16a)]. Since the crystal may 
lose heat by natural convection and/or thermal radiation, one has to 
resort to theoretical relations in fluid dynamics 29 and radiative trans- 
fer, 18 respectively, for guidance in the selection of h. For the purposes 
of this paper, the coefficient for heat transfer by a gas from a vertical 



Table I — A summary of input parameters 



2 


4 


0.2 


7.2 


0.0005 


0.003 


200 


7200 



Lower Upper 

Property Limit Limit Fixed 

Radius, ro(cm) 

Length, / = pt(cm) 

Pull rate, p(cm/s) 

Time, t(s) 

Melting point, 7>(°C) 1238 

Average ambient temperature 100 300 

drop, T, - r„(°C) 
Thermal conductivity, 0.08* 

Jf{watts/cm K) 
Thermal diffusivity, ic(cmVs) 0.04 

Thermal expansion coefficient, 1 x 10~ 5 

a(°K-') 
Elastic stiffnesses, C, i (dynes/cm ) 10 x 10 " 

C l2 (dynes/cm 2 ) 4.6x10" 

C« (dynes/cm 2 ) 5.1 x 10" 

* Values in column are near Tf. 
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wall can be written in the condensed form 

1/4 



*[cm- l ] = A'^^) P l/2 , 



(42) 



where I is the length of the wall, P is the total pressure, and h is a 
function of the thermal and transport properties of the gaseous am- 
bient and the solid wall. The expression for a liquid medium is similar 
to eq. (42) with the single exception of taking P = 1. A critical 
evaluation of the constants required to estimate h! will be given 
elsewhere. 28 

Based on the Stefan-Boltzmann equation for thermal radiation, the 
coefficient for heat transfer by radiation can be obtained from 16, 18 

Mem- 1 ] = 5.672 X 1(T 12 e ~ Ta ) X (43) 

where c and Jf are the total emittance and thermal conductivity of 
GaAs, respectively. We have recently evaluated e of n-type GaAs as a 
function of temperature, doping level, and thickness by an analysis of 
the absorption spectra. For a crystal of 2-cm radius containing 10 16 
cm -3 n-type impurities, e increases from 0.52 to 0.57 in the temperature 
range between 1200° and 1400°K. 30 If (T - T a ) is not very large, the 
temperature-dependent term in eq (43) can be given within a good 
approximation by AT'i to yield for h i6 

Mcm-'] = 0.0227( T5 ^)W (44) 

In Table II, the convective and radiative heat transfer coefficients 
are summarized at three ambient temperatures. To perform the cal- 
culations according to eqs. (42), (43), and (44), we assumed that (7>- 
T)/(Tf — T a ) = 0.7 and 0.4 for gaseous and liquid or radiative heat 
transfer, respectively. We can immediately conclude from Table II 
that in the case of lec growth for T, - T a > 100, natural convection 
via the B 2 3 (Z) competes with radiative transfer in dissipating the heat 
from the crystal's external boundaries. This unexpected result is due 
to the relatively low viscosity of B 2 3 at these temperatures (100 poise 
at 1000°C). 31 Hence, convective transfer via the gas phase is nearly 
negligible in comparison with that by the encapsulating liquid. Fur- 
thermore, we find that the approximate eq. (44) for h yields the value 
of the radiative heat transfer coefficient with adequate precision. Since 
the convective and radiative heat transfer coefficients are additive, h 
a 0.6 is an appropriate approximate value for lec growth if T f — T a 
= 200°K. 
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4.2 The effect of the variables on the dislocation density 
4.2. 1 Location and heat transfer coefficient 

The effects of axial position and h on the o,„, pattern are illustrated 
in Figs. 7 and 8, respectively. In view of the fourfold symmetry of the 
pattern demonstrated in Fig. 6, a polar plot may represent one axial 
location per quadrant. In Fig. 7, four positions between the top of the 
crystal and 0.2 cm from the solid-liquid interface are shown, assuming 
that h = 0.6 cm"', r = 2 cm, andp = 0.001 cm/s and that 5400 seconds 
have elapsed since the beginning of growth. The numerical labels on 
all the contours translate immediately into absolute a to t values if 
multiplied by (aE/1 - v)[(T f - T«)/200]. Since this convention is 
observed throughout the paper, the effect of a change in ambient 
temperature is essentially reflected by a change in 7/ — T a . 

One can readily conclude from Fig. 7 that at a given time the stress 
rises from a low value at the top to a maximum at about 0.6 cm from 
the solid-liquid interface and then diminishes as the interface is ap- 



1.6 cm FROM 

TOP 

5400 sec 




4.8 cm FROM 

TOP 

5400 sec 



5.2 cm FROM 

TOP 

5400 sec 



Fig. 7— The axial location dependence of o,„, or dislocation density contour lines for 
an <001> GaAs boule at t = 5400 s and h = 0.3 cm '. The o to t levels are labeled, which 
can be converted into absolute stresses when multiplied by (aE/1 - r)[(Tf — T„)/200]. 
Each successive quadrant in a counterclockwise direction illustrates the dislocation 
distribution in a wafer at a progressively shorter distance from the solid-liquid interface. 
The relevant parameters are p = 0.001 cm/s, length = 5.4 cm, r = 2 cm, and k = 0.04 
cm 2 /s. 
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Fig. 8— The dependence of o„„ or dislocation density on the heat transfer coefficient 
(h) and location (-z) for a <001> GaAs crystal at the center and <100) edge of a wafer. 
Absolute stress values can be obtained by multiplying the ordinate with (aE/l — i»)[(7/ 
- T a ) /200]. The labels correspond to the distance in centimeters from the top of the 
crystal. Other key parameters are t = 5400 s, length = 5.4 cm, r = 2 cm, and k = 0.04 
cm /s. The curves for the center at the top and 1.6 cm from the top, which should be 
shown merging with the 3.2-cm line at low h values, were terminated at h = 0.06 cm" 1 to 
reduce crowding. 

proached. Moreover, it can be seen that the most severe stress gradient 
is at the periphery. 

In Fig. 8 we show the dependence of a to t on h at the center as well 
as at the (100) edge (r = 2 cm, r = 2 cm) of a {100} wafer. Except for 
the additional slice at 3.2 cm from the top, the other four locations are 
the same as in Fig. 7. On the whole, a tot monotonically rises with h at 
both the center and edge of the wafer. However, locations in the upper 
part of the ingot (seed-end) exhibit a maximum at an h value which 
increases with increasing distance from the top. The simplest mathe- 
matical explanation of this result follows from the properties of the 
term e'* 3 " 4 ', which plays a prominent role in the limit of large \p, and 
moderate /?„ in the expressions for T — T a /T f — T a [eq. (28)] and a tot 
[eq. (39) via eq. (32)]. As h becomes larger, so does /?„; hence, for any 
\p eventually e~ fi "* -* 0. The smaller the value of \p (the closer the cut 
to the interface), the bigger the h at which the downturn occurs. 

A more physical interpretation of the maximum in the a to i versus h 
curves is as follows: The crystal growth model we have been investi- 
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gating includes two extreme cases. On the one hand, when h is very 
small, T-T a /Tf-Ta-*l, because the crystal temperature uniformly 
approaches that of the molten GaAs(7». On the other hand, when h 
is extremely large T - T a /T f - T a -+ because the ingot reaches the 
ambient temperature (T a ) almost instantly. Regardless of the limiting 
case [(T - T a )/(T f - T a ) = or 1], the thermal stresses must diminish 
if the crystal temperature is uniform and thus possesses a maximum at 
an intermediate value of h. Moreover, the temperature at the top of 
the crystal is closer to T a than that of the middle. This explains for the 
same h the tendency of the stress at the top to decrease while at the 
middle it is still increasing. 

Similar to Fig. 7, Fig. 8 also shows the maximum a to t at 0.6 cm from 
the solid-liquid interface for both the center and edge of the wafer. At 
this location, a to t is a monotonically increasing function of h in a range 
between 0.01 and 1. Thus the maximum dislocation density is attained 
at a relatively short distance (say, <1 cm) from the interface. Conse- 
quently, by lowering h near the interface, a reduction in dislocation 
density can be accomplished. 

4.2.2 Radius 

In Fig. 9, the a to t contour lines for a crystal with 2-cm radius are 
superimposed over the lines for one with a 4-cm radius. The other 
parameters are h = 0.6 cm -1 , t = 4500 s, p = 0.001 cm/s, and four 
locations between the top and 0.5 cm from the interface are given. In 
general, considering identical locations, the edge of the crystal with a 
larger radius is more highly stressed than the one with a smaller 
radius. Comparing equivalent points with the same p = r/r , the effect 
of increasing the crystal diameter on o to t appears to be superlinear. 
However, for cuts near the interface (0.5 cm) the stresses at equivalent 
points are almost the same. This is a consequence of the fact that, at 
»// = 0.5 cm, the a to t of a small diameter crystal just reached its 
maximum, whereas that of a large diameter crystal is already over its 
maximum. 

To a great extent, the variation of ow with r mimics the variation 
with h (Fig. 8) because the dimensionless quantity hi = r h occurs in 
the key equations. But r also appears implicitly in ^ = (z + pt)/r. 
Hence, any h curve family such as Fig. 8 can be used to visualize the 
effect of increasing the radius provided that appropriate adjustments 
are made in the hi and \p coordinates. A key finding which is consistent 
with Fig. 9 is that the larger the diameter of the crystal the farther 
from the interface is its location of maximum stress. 

Although crystals with a larger r are exposed to higher stresses, 
there is an inherent advantage in growing such crystals if discarding a 
portion is an acceptable production practice. For example, in the third 
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Fig. 9 — The dependence of o lol or dislocation density contours on radius and axial 
location at 4500 s and h = 0.6 cm" 1 . The inner and outer circles confine the profiles for 
r = 2 and 4 cm, respectively. The labels and labels with bars denote identical stress 
levels in the GaAs crystal of 2- and 4-cm radius, respectively. Multiplication of the levels 
by (aE/l — v)[(Tf— T„)/200] yield the absolute stresses. Each successive quadrant in a 
counterclockwise direction depicts the dislocation distribution in 2- and 4-cm radius 
wafers at a progressively shorter distance from the solid-liquid interface. Additional 
parameters are p = 0.001 cm/s, length = 4.5 cm, and k = 0.04 cm 2 /s. 

quadrant of Fig. 9 (3.2 cm from the top), a sizable quasirectangular 
area of the 4-cm radius crystal is at a to t ^ 56. In the case of the 2-cm 
radius crystal, almost the entire area is at a comparable level or less. 
If level 56 were below the threshold stress for slip, then there would be 
an approximately 2:1 gain in area yield in cutting out a rectangular 
portion of the larger diameter wafer over keeping the entire area of 
the smaller diameter one. Since we have found that the stresses in the 
slender ingot would further increase down to a distance of =*0.5 cm 
from the interface (fourth quadrant), the choice to cut a sizable area, 
low in defect density, from the broader ingot is even more favorable. 

4.2.3 Time 

Time is an implicit variable in the quasi-steady-state temperature 
and stress profile equations. In the preceding parts of this section, the 
stress levels were given at various axial locations, subsequent to the 
growth of a GaAs ingot for a fixed length of time. In Fig. 10, we present 
polar plots of the a to t contours which prevail at the top of the crystal — 



DISLOCATION GENERATION IN GaAs CRYSTALS 621 



TOP. 1.8 cm 
LONG 

1800 sec 



TOP, 4.5 cm 
LONG 

4500 sec 




*- <100> 



0.5 cm 
LONG 
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Fig. 10— The time dependence of the a,„, or dislocation density contours of a GaAs 
boule at h = 0.6 cm -1 . Starting in the fourth quadrant and proceeding in a clockwise 
direction, the time evolution of the dislocation pattern at the top end of a growing 
crystal is shown. The labels correspond to the stress levels and are converted into 
absolute values when multiplied by (aE/1 - v)/[(T f - 7 1 „)/200]. Other relevant param- 
eters are p = 0.001 cm/s and r = 2 cm. 

a fixed location— after 500, 900, 1800, and 4500 seconds of growth. The 
other key parameters are h = 0.6 cm" 1 , r = 2 cm, andp = 0.001 cm/s. 
One can readily see that, as crystal growth proceeds, the top is exposed 
very early to the maximum stress, namely between 500 and 900 
seconds.* Afterwards, the stress levels rapidly decline in magnitude. It 
is apparent that the patterns at a fixed location but for different 
durations of growth (Fig. 10) parallel the ones at a set time but for 
different axial locations (Fig. 7, Fig. 9, r = 2 cm). Intuitively, this 
should not be surprising. For instance, after 4500 s, the top is 4.5 cm 
from the solid-liquid interface and a segment of the crystal is at a 
distance of, say, 0.5 cm from the interface. Then the pattern at the 
latter position may resemble that of the top at 500 s, at which time the 
boule was only 0.5 cm long. Formally, the time-location equivalence is 
only exact when e _/w is a very good approximation of the axial function 
in eq. (28). Otherwise, it indicates a useful trend in correlating time 
and position variations. 



* Compare Fig. 10, quadrants 3 (900 s) and 4 (500 s) near the center. 
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By concentrating on a single point at the top, a more complete 
representation of ct (0 i can be provided than has been achieved by the 
polar plots alone in Fig. 10. We show in Fig. 11 a to t as a function of 
time at the top (100) edge of a GaAs ingot. A family of these curves 
with h varying between 0.01 and 2 cm" 1 is given. Each member of the 
family exhibits a maximum in a to t at a critical time, t c . More effective 
heat transfer in the growth system leads to an earlier t c and to a larger 
and steeper maximum in a tot . As already mentioned with respect to 




3000 4000 5000 

TIME IN SECONDS 



6000 



8000 



Fig. 11— The variation of o tu , or dislocation density with time and heat transfer 
coefficient (h) at the top ( 100) edge of a GaAs crystal. The labels are the heat transfer 
coefficients in cm '. Multiplication of the stress levels by (aE/\ - v)[(T/ - T„)/200] 
yields the absolute value of the stresses. The other parameters are p = 0.001 cm/s, r 
= 2 cm, and k = 0.04 cm"/s. 
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Fig. 8, it can also be seen in Fig. 10 that as h -> and oo, a tol -> 
because the temperature of the crystal approaches uniformly the 
constant values 7} and T a , respectively. 

If h = 0.6 cm -1 , the maximum in a to t at the top <100) edge of GaAs 
is at ~700 s. As the crystal grows longer, the size of the stress 
diminishes. Thus, the time-dependence study confirms the conclusion 
that the maximum dislocation density is attained during the early 
phase of growth. Of course, it is presupposed here that the dislocations 
form instantaneously when the stress is above the critical resolved 
shear stress and that the density of dislocations thus frozen-in irre- 
versibly will not be reduced by the drop in the stress level as the 
crystal grows longer. A downward adjustment of h is the most obvious 
method by means of which the dislocation density is reducible. 

4.2.4 Pull rate 

Figure 12 presents a composite set of polar plots to demonstrate the 
effect of pull rate on a, ot . The left and right halves show the results for 



<110> 



TOP, 1 
LONG 
3600 sec 
600 sec 




TOP. 1.8 cm 
LONG 

3600 sec (U 
600 sec (L) 



Fig. 12— The effect of pull rate (p) on the o,„, or dislocation density contours at the 
top of a GaAs boule. The labels and labels with bars (space permitting) denote identical 
stress levels at p = 0.0005 and 0.003 cm/s, respectively. Multiplication of the levels by 
(aE/l - v)[(T f - T„)/200] provides the absolute stresses. The two quadrants in the 
upper and lower halves show the patterns at h = 0.6 and 0.06 cm ', respectively. The 
two quadrants in the left and right halves illustrate the distribution at the top of a 1.8- 
and 3.6-cm long crystal, respectively. Since the lengths are constant, the times for 
growth are consistent with the selected pull rates. Other parameters of interest are r 
= 2 cm and k = 0.04 cm 2 /s. 
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the top of a 1.8- and 3.6-cm long ingot, respectively, while the top and 
bottom halves represent the profiles at h = 0.6 and 0.06 cm -1 , respec- 
tively. The other parameters are r„ = 2 cm and p = 5x 10 -4 and 3 x 
10~ 3 cm/s. The fixed lengths correspond to growth times t = 3600 s 
and 600 s forp = 5 X 10" 4 and 3 X 10~ 3 cm/s, respectively, for the 1.8- 
cm long ingot. When the length is 3.6 cm, the times are doubled. As 
one would expect from the analytical expressions and as also seen in 
Fig. 12, the influence of a sixfold increase of pull rate on the dislocation 
distribution is more pronounced in the case of a longer boule than that 
of a shorter one. 

To provide a more precise discussion with reference to the effect of 
pull rate on a tot , an approximate expression will be derived for the 
fractional increase of o, 0( at the top (100) edge as p changes. The 
calculation can be performed solely by means of eq. (32b) since at the 
(100) edge a iot = (8w/6)<r#. Then, in view of the fact that pi is more 
effective in modifying h p than /?„,* one finds for the relative change in 
ow, (otot — atot)/<Ttot-, taking only the leading term in the series of eq. 
(32b), 

„" ~' 

Otot Otot 



O tot 



-(pr-p'i)^- 



1 _ p-WAo 

1- 



l/r [h P ' + /?. - (K - 0i)e- w/r °] 



(45) 



where the double and single primes refer to faster and slower pull 
rates, respectively, and I is the length of the crystal. At a reasonably 
rapid h, eq. (45) reduces to 

Otot — o'tot ... , . / / 1 \ lMM 

=* (Pi -PihH 1 — ; (46) 



°' t01 2r °\ l ih" M /? \/ 

since the exponentials approach zero. For very small h\ , a Taylor series 
expansion of the exponentials yields 

- (Pi -Pih-l 1 ; • (47) 



ro 

According to eq. (45), the effect of stress can be readily observed for 

a long crystal if the difference in pull rates is large. If the heat transfer 

coefficient is substantial, the full impact of the growth rate term 

e Pt+/2 m e q (32b) i s achieved because the multiplier 1 - [l/l/r (hp + 



* Eas. (29) show that (i„ is relatively insensitive to a change in p, on account of its 
dependence onp'f/2. 
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/*i )] -* 1 [eq. (46)]. On the other hand, if h is minute the multiplier 
approaches zero [eq. (47)] and the dependence of a to t on pull rate, 
especially for short I, is negligible. The stress patterns in the top two 
quadrants (h = 0.6) and bottom two quadrants (h = 0.06) of Fig. 12 
fully demonstrate the validity of these general remarks. 

It is apparent that, of all the parameters here considered, a change 
in pull rate is perhaps the least effective method by which one may 
influence the stress levels in GaAs. As can be most readily observed in 
the first quadrant of Fig. 12, for identical spots the wafer sliced from 
an ingot grown at a slower pull rate is under less stress than the one 
obtained at a faster pull rate. This can be rationalized by noting that 
the key quantity h p is the sum of h\ andpi/2. In essence, a more rapid 
pi is analogous to an increase in hi ; hence, the greater pi , the steeper 
the temperature gradient from axis to edge and the larger the thermal 
stresses. Although the effect of pull rate on the stress pattern is 
relatively small, if on account of other factors the resolved shear stress 
happens to be only marginally higher than its critical value, it may 
become feasible to reduce the dislocation density by lowering the pull 
rate. 

4.3 Comparison with experiment 

In Fig. 13 a macrophotograph of a {100} wafer cut from the vicinity 
of the seed or top-end of a Te-doped GaAs boule is presented. To 
reveal the dislocation structure, the wafer was etched in fused KOH at 
300°C for 1 hour. It has been previously reported by Angilello et al 32 
that, on a {100} wafer of GaAs, there is a one-to-one correspondence 
between the dislocation density determined by X-ray transmission 
topography and the pits decorated by the etchant where the disloca- 
tions intersect the surface. The dislocation density in Fig. 13 varies 
between 10 4 and 5 X 10 4 /cm 2 . 

A microscopic view of the dislocation distribution for a Cr-doped 
{100} wafer is given in Fig. 14. This composite picture clearly illustrates 
the dislocation density in four critical regions of an etched GaAs slice. 
One can readily discern the following key features of the dislocation 
distribution: 

(i) Fourfold symmetry. 

(ii) Maximum density at the (100) edge. 

(Hi) Minimum density midway between center and (110) edge. 

(iv) Intermediate densities at center and (110) edge, but the edge 
density is somewhat higher. 

In fact, all the calculated o to t contour lines possess the preceding 
characteristics. To make the comparison more convenient, we have 
replotted the polar diagram of Fig. 6 in the form of a "gray scale" in 
Fig. 15. Accordingly, increasing stress levels in Fig. 6 correspond to 
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Fig. 13 — Macrophotograph at ~5X magnification of a KOH-etched {001} GaAs wafer. 
The crystal was cut from the vicinity of the top-end of a Te-doped GaAs boule grown by 
the lec technique. 



more darkly shaded areas in Fig. 15, in the manner of a geographic 
contour map. The agreement between the experimental results and 
the theoretical prediction is very good. 

The observed symmetric pattern in Figs. 13 and 14 is not dependent 
on the type of dopant incorporated in GaAs in view of the similar 
distribution of dislocations found in chromium as well as tellurium- 
doped samples. However, the farther away the cut from the top of the 
ingot the more diffuse the distribution. This is not surprising if account 
is taken of the tendency of dislocations to move out of their slip planes 
by climb. Morever, the typical 60-degree dislocations 33 generated at 
the top by thermal stress-induced glide continue into the next-to-grow 
layer of the crystal and add to the glide dislocations arising therein. 
Therefore, extra mechanisms notwithstanding, the progenitor of dis- 
locations in GaAs is crystallographic glide relieving the severe thermal 
stresses associated with the Czochralski growth process. 
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Fig. 14 — Photomicrographs of important regions within a KOH-etched {001} GaAs 
wafer. The crystal was cut from the vicinity of the top-end of a Cr-doped GaAs ingot 
grown by the lec technique. The original locations are indicated in a realistic though 
schematic fashion. 



4.4 Estimation of the critical stress level 

Having correlated the empirical dislocation density pattern of a 
{100} GaAs wafer with the total resolved shear stress profile, a to t, the 
critical value of this quantity at which slip is initiated becomes a topic 
of preponderant interest. It has been previously stated that in all 
diagrams multiplication of the stress level by (aE/1 — v)[(Tf — T a )/ 
200] yields the stress in absolute units (dynes/cm 2 ). Although the 
original stress equations are those of isotropic thermoelasticity, by a 
minor modification of the constant E/l — v it is possible to take into 
account to some degree the anisotropy of the cubic GaAs. According 
to Brantley, for stresses acting in directions within {111} planes E/l 
- v is an invariant quantity. 34 Since the GaAs slip system is {111}, 
( 110), it is reasonable to employ (E/l —v) {mj in our estimates. Based 
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Fig. 15 — "Gray scale" representation of the constant otot or dislocation density contour 
lines for the top wafer of a (001 > GaAs boule. The shading is based on the pattern in 
Fig. 6. 



on Brantley's work we can write 
E 



i - v) S? n + Sfa - £f/3 

{in} 



or 



E 



1 - v 



= C, + Ci2 + 



(111) 



C(C n + 4C 12 ) - 3C? 2 
2C + 1.5C n 



(48a) 



(48b) 



depending on whether one uses elastic stiffnesses (Cy) or compliances 
(Sfij). The symbols ,^and C represent the departures from isotropy and 
are given by Sf= ^ n - ^12 - KSfu, C = C u ~ %(C„ - C12). If S/>= C 
= 0, the crystal is isotropic. 

Table III lists the expansion coefficient, (aE/1 — p){ih), (aE/1 — 
»*)iBotropic for GaAs at various temperatures. The data sources and their 
critical evaluation will be given elsewhere. 28 Obviously, the product of 
(Tf— T a ) /200 times the tabulated value times the stress level given in 
the various figures provides the absolute stress. 

The critical resolved shear stress (crss) of GaAs between 250° and 
550 °C has been determined by Swaminathan and Copley. 35 To extrap- 
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Table III — Some elastic properties of GaAs 



Temperature (°K) 


T, 


7>- 100° K T, 


- 200°K T, 


-300°K 


298°K 


a x 10- 6 ( o KT') 


10.4 


10.1 


9.69 


9.31 


5.80 


( aE \ k 












U-'J C1I1) 

10 7 (dynes/cm 2o K) 


1.55 


1.51 


1.48 


1.44 


1.01 


l^E- ) x 10 7 (dynes/cm 2o K) 


1.10 


1.07 


1.05 


1.02 


0.72 


( aE \ 

10 7 (dynes/cm 2o K) [aniso- 
tropic coefficient] 


1.32 


1.28 


1.26 


1.22 


0.86 


8 x crss x 10 7 (dynes/cm 2 ), 


4.5 


5.2 


6.2 


7.5 





olate their results to elevated temperatures, we assumed, as is usual, 25 
that the crss is an exponential function of reciprocal temperature. 
Table III also lists crss multiplied by 8, which we call the effective 
crss. The factor 8 is used to facilitate comparison with a to t because at 
the (100) edge of a {100} wafer 8 slips of identical magnitude operate 
[atot = 8v / 6/6a < ,]. 

Let us compare Fig. 11 with the values in the T f - 200°C column in 
Table III. It can be seen that even for level 5 o to t = 5 X 1.5 X 10 7 = 7.5 
X 10 7 dynes/cm 2 is larger than the effective crss of 6.2 X 10 7 dynes/ 
cm 2 . Therefore, if the plotted level is above 4.1, the (100) edge of the 
top will develop dislocations at -1000 s. This corresponds to a thresh- 
old of h = 0.02 in Fig. 11. In other words, if h< 0.02 cm -1 dislocations 
cannot be avoided in GaAs at least at the (100) edge of a 4-cm 
diameter {100} wafer. At h = 0.6 cm" 1 which is the heat transfer 
coefficient appropriate for lec growth, the estimated atot variation in 
the (110) direction is -21, 12, and 54 X 10 7 dynes/cm 2 from the center 
through a midway point, and the periphery, respectively (Fig. 10, 
fourth quadrant). At the (100) edge ow = 105 X 10 7 dynes/cm 2 (Fig. 
11). Although all these values are above the effective crss, one should 
be aware of the errors involved in the various estimates and realize 
that it may be relatively easy to achieve a nearly dislocation-free 
material near the (110) midway spot (12 vs 105 X 10 7 dynes/cm 2 ). 
Therefore, by judicious design changes to decrease Tf — T a , h and p 
which are associated with the Czochralski growth process of GaAs, at 
least the reduction, if not the elimination, of dislocations is a realistic 
possibility. Since the dislocation density gradient is the steepest near 
the periphery, even improvements in the inner two-thirds of the wafer 
area could be of practical benefit. 
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4.5 The effect of anisotropy 

To calculate a to t, we adjusted the isotropic thermal stress equation 
by using an invariant value of E/l — v valid for stress directions within 
{111} planes. 34 Strictly speaking, this is not identical with finding a 
true anisotropic solution for the thermal stresses arising during Czo- 
chralski growth. Even under isothermal conditions, to solve the equa- 
tions of elasticity for an anisotropic body is a formidable undertaking. 36 
Only in the case of a simple temperature profile is there any hope to 
obtain in a closed form the thermal stresses acting on an anisotropic 
cylinder. This can be most clearly seen by inspecting the partial 
differential equation of thermoelasticity for a cubic crystal with a 
(100) axis, in the so-called stress formulation. Following the derivation 
outlined by Boley and Weiner 37 for plane strain in an isotropic body, 
we find for cubic material 



SfwSf (d 2 o x d 2 o y \ dS*T 

<? 2 u-9> 2 Adx 2 + dy 2 ~ Sh + Sk 



V'(a, + o y )- — 2 -A — , + — i = - , - , (49) 



where V 2 is the Laplacean operator in Cartesian coordinates. If £f= 0, 
the partial differential equation is that of an isotropic body and the 
complication due to the inclusion of the second term on the left-hand 
side of eq. (49) disappears. Grechushnikov and Brodovskii 38 succeeded 
in solving eq. (49) for a quenched cubic cylinder with a (100) axis by 
assuming a radially parabolic temperature profile taken from the early 
heat transfer studies of Adams and Williamson. 39 

Based on the work of Grechushnikov and Brodovskii, 38 we can 
evaluate the ratio (^(anisotropic/anisotropic) in the form 



ag( anisotropic) 



(50) 



^(isotropic) 1 - [^ii^74(5i, - y 2 2 )] 

and substituting the appropriate values 

^(anisotropic) = 1.2a e (isotropic). 

Since at the (100) edge of a crystal ff, ot = 8(>/6/6)<70, multiplication of 
the values in the values in the third column in Table III [(aE/1 — 
yjiso] by 1.20 could serve as a suitable estimate of the effect of 
anisotropy on a to t . As shown in Table III, the anisotropic coefficient 
lies between the (aE/1 - w)[m) and (aE/1 — i>) iso results. In practical 
terms, this means that the threshold level rises from 4.1 to 4.9 if the 
anisotropic coefficient instead of (aE/1 — J»){in) is used. This is indeed 
a minor perturbation and the conclusions of this study are not influ- 
enced by the anisotropic solution. 
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4.6 Additional characteristics of the quasi-steady-state solution 

The temperature profiles prevailing during Czochralski growth [eq. 
(28)] were based on the quasi-steady differential equation [eq. (5)] 
which depends explicitly on pull rate and implicitly on time, via \p and 
\p, [eqs. (29)]. Since T a was taken as a constant, in the limit of p -> 0, 
the qss solution reduces to the steady-state result of Brice 14 /Carslaw 
and Jaeger. 16 To gain insight into the degree of approximation involved 
in using the QSS instead of a true time-dependent solution of eq. (4) 
when determining the temperature distribution, we have examined the 
departure from equality of eq. (4) by substituting T from eq. (28). 
Having performed the indicated operations, we conclude that at the 
top-end of the crystal the qss solution satisfies eq. (4) if and only if 
exp(pi/2 - p n )pt/r approaches zero. Since, in general, pi/2 < /?„, this 
finding suggests that a moderate pull rate, efficient heat transfer into 
the ambient medium, and a reasonable passage of time since the 
inception of growth are the conditions at which the qss solution is the 
most satisfactory. 

Not only mathematical but also conceptual difficulties are associated 
with a search for a true transient solution to Czochralski growth. In 
classical heat transfer (or diffusion) problems, a physical body exists 
initially (t = 0) with given properties (prescribed temperatures or 
fluxes). In contrast, in crystal growth, the initial state is that of a liquid 
from which the solid materializes by phase transformation at t > 0. 
Hence, the boundary conditions for the crystal [e.g., eqs. (15) and (16)] 
have no initial relevance. 

Strictly speaking, heat transfer during crystal growth belongs to the 
group of problems involving a change of state which is often referred 
to as "Stefan's problem." According to Carslaw and Jaeger, 16 the few 
existing exact solutions of "Stefan's problem" pertain to very simple 
boundary conditions and semi-infinite geometries. A typical illustrative 
example is as follows: Suppose a liquid column at T f is contained in a 
very wide and tall crucible. Then, beginning at t > the temperature 
at the bottom is lowered to a fixed value below 7> and, consequently, 
a freezing front moves upwards. Clearly, each region must satisfy the 
one-dimensional partial differential equation of heat conduction and 
the temperatures are equal at the solid-liquid interface (T f ). In addi- 
tion, if the density of the two phases is nearly equal, the heat-flux 
balance at the interface demands that 

x «. < m + ^:* (Bi) 

dS dS V m dt 

where s is the distance from the bottom and Xi(i = crystal or liquid), 
bH f , and V m denote the thermal conductivity, heat of fusion, and 
molar volume, respectively. In spite of the nonlinear boundary condi- 
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tion [eq. (51)], Neumann succeeded in determining the exact vertical 
temperature profile in both phases and the location of the interface as 
a function of time in terms of error functions. 16 

Currently, more sophisticated situations such as a body of finite 
length can be treated by Boley's "embedding technique" which in- 
volves an imaginary extension of the solid and liquid phases into the 
space actually occupied by the diminishing liquid and growing solid 
fronts, respectively. 40 Although this technique circumvents the conflict 
between the boundary and initial conditions, it leads to a complicated 
integro-differential equation which is only amenable to a series solu- 
tion. 

Therefore, it is apparent that only by overcoming the severe con- 
ceptual and analytical obstacles can one attain a more realistic descrip- 
tion of the Czochralski growth process especially in the vicinity of the 
solid-liquid interface and at the beginning of crystallization. In partic- 
ular, if the effect of the physical and geometrical conditions on the 
interface shape is a matter of preponderant interest, a numerical finite- 
difference solution of the appropriate heat transfer equation 41 is likely 
to be unavoidable. 

As a reasonable approximation for GaAs under the prevailing Czo- 
chralski growth conditions, it has been assumed in the qss model that 
at the interface (\p = s/r a = 0) T = Tf. Below Tf in the crystal the qss 
isotherms [\j/ = f(p, T = constant)] are all concave with a maximum at 
p = 0. Essentially, the isotherms mimic the complement of the tem- 
perature profiles in Fig. 3 [1 - (7> - T)/(7> - T a )]. Of course, at the 
interface a flux balance of the type given in eq. (51) must be obeyed. 
In general, depending on whether the temperature of the liquid nutri- 
ent is higher or lower than on the axis at Tf, a convex or concave 
interface, respectively, results. 42 However, it can be demonstrated with 
respect to the qss calculations here considered that the assumption of 
a planar isotherm at Tf is internally consistent. Let us evaluate the 
heat flux according to eq. (51). Then, taking dTi/ds = and identifying 
ds/dt = p, 

fiT Aff 

X-r 1 =p-rr 1 = -3.7 watts/cm 2 . (52) 

ds V m 

Differentiating eq. (28) and substituting the appropriate parametric 
values,* we find the fluxes at 0.4 cm from the solid-liquid interface and 
at the interface at selected values of h as shown in Table IV. The listed 
results clearly indicate that if h = 0.1 or less, nearly planar interface 
can be maintained; moreover, the heat of crystallization is readily 



* bHf = 25.080 kcal/gmole,"" p = 5.16 g/cm" 4 -", / = 3600 s, length = 3.6 cm, p = 0.001 
cm/8, r = 2 cm, .rf= 0.08 watts/cm°K, k = 0.04 cnr/s, T,-T„ = 200°K. 
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Table IV 

/i[cnT'] jr c — ^[watts/cm 2 ] 

OS 



3.2 cm from top 3.6 cm from top 

axis/edge axis/edge (interface) 

0.03 -1.8/-2.2 -1.8/-3.4 

0.10 -4.0/-5.1 -4.1/-9.2 

0.30 -6.7/-8.9 -7.1/-20.1 



dissipated if h S 0.1. It should be emphasized that irrespective of the 
shape of the melting point isotherm, its effect on the temperature 
distribution at some distance from the interface — where the maximum 
in atot and the dislocation density occurs— is not expected to be 

significant. 

An additional feature of the qss model is that the more effective the 
convective heat transfer, the more concave the isotherm. The partial 
derivatives of T obey the relationship 

'ds\ -{dT/dp)s 



— = 



dp J (dT/ds), 



(53) 



at any isotherm. Evaluating eq. (53) by means of information such as 
given in Fig. 3 and Table IV, we find at 0.4 cm from the interface and 
at p = 1 



/i[cm-'] 



0.1 
0.3 



%)r 



-6.9 
-9.3 



Obviously, the more negative tangent at h = 0.3 corresponds to a more 
concave isotherm. 

V. SUMMARY AND RECOMMENDATIONS 

We have determined the temperature distribution arising during the 
Czochralski growth of GaAs by solving the quasi-steady-state partial 
differential equation for heat conduction. The solution includes time, 
radius, axial location, pull rate, a constant ambient temperature, and 
heat transfer coefficient among the variables and permits the assess- 
ment of their effect. The temperature profiles enabled us to calculate 
the radial, tangential, and axial thermal stresses in the growing cylin- 
drical boule in a closed form. These stresses were required to evaluate 
the sum of the absolute values of the 12_resolved shear stress compo- 
nents (atot) appropriate for the {111}, (110) slip system. It was postu- 
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lated that in the case of complete stress relief o to t is proportional to the 
dislocation density within an additive constant. To directly relate the 
model calculations to actual sample geometry, we have generated 
constant o,„, or dislocation density contour lines for a circular {100} 
wafer. The dislocation pattern has fourfold symmetry. An absolute 
minimum in density is found in the (110) direction at ~0.6x wafer 
radius. The center and edge of the wafer are heavily dislocated, the 
density being largest at the edge. The calculated patterns are in accord 
with KOH-etched dislocation patterns on wafers cut from near the top 
end of Cr- and Te-doped GaAs ingots, strongly suggesting that the 
original source of dislocations is crystallographic glide induced by the 
excessive thermal stresses associated with crystal growth. 

By a variation of the parameters, we have reached the following 
conclusions: 

(i) Doubling the radius (2 to 4 cm) of the ingot more than doubles 
the dislocation density. However, for comparable densities, the larger 
diameter wafer has a useful area in excess of the entire area of the 
smaller diameter one. 

(ii) Reducing the pull rate by a factor of 6 (10.8 cm/h to 1.8 cm/h) 
has only a small effect on reducing the dislocation density. 

(ill) Dislocations near the top of the ingot are grown-in by the time 
the ingot is less than 1 cm long (i.e., a lo t is at its time-dependent 
maximum). 

(iv) The more effective the convective and radiative heat transfer 
from the surface of the boule, the larger a to t and the dislocation density. 
In the case of lec growth, natural convection (h = 0.24 cm -1 ) and 
radiation {h = 0.34 cm" 1 ) through the ~l-cm long B2O3 column are the 
predominant heat transfer mechanisms. For Czochralski growth with 
a gaseous ambient, h is only about 2 to 4 percent of that for B2O3. 

(v) Using extrapolated values of the measured elastic constants, 
thermal expansion coefficient and critical resolved shear stress in 
combination with the model calculations, we find that dislocations can 
be avoided if the approximate inequality (T f — T a ) h?s 4°K/cm holds. 
For (T f - T a ) = 200°K, h should be less than 0.02 cm" 1 . Therefore, it 
is not surprising that the lec growth of GaAs (h a 0.6) leads to 
significant dislocation generation. Although, in practice, h cannot be 
readily decreased below 0.02 cm -1 , any reduction in h and/or Tf — T a 
will result in a proportional reduction in dislocation density. 

Some of these conclusions appear to be consistent with the experi- 
ence of Leung and Allred 40 who reported on a liquid seal growth 
technique for GaAs. In this method, B2O3 is located at the growth 
chamber pull-rod junction. The chamber is filled with excess AS2 and 
a resistance-wound after-heater over the ingot is employed. That this 
technique provides, as reported, low-dislocation density crystals is 
reasonable in view of the quasi-steady-state model. Accordingly, based 
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on molecular weight considerations h for As 2 is very small (must be 
below that for A in Table II), convection through B 2 3 is avoided, and 
(T f - T a ) is reduced by means of the after-heater. However, for routine 
growth of GaAs, the above method is impractical. 

Clearly, there are several areas to focus on in the case of lec growth. 
The thinnest possible B 2 3 layer consistent with preventing As escape 
should help limit the large connective transfer well before the dislo- 
cations reach their maximum density. Radiation shields will lower the 
radiative heat transfer coefficient and after-heaters may be used to 
increase T a . Moreover, on large area wafers, the high dislocation 
density periphery may be removed following growth. In any event, an 
empirical balance must always be struck between reducing the dislo- 
cation density by making the system more isothermal and the antici- 
pated difficulty to pull usable crystals in a low temperature gradient 
environment. 
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